#!/usr/bin/python
# GM_darklight.py
'''Compares the distribution of (potentially unseen) compact matter, from Kroupa, to visible matter, from Binney & Merrifield.

It is sufficient to compare the two distributions locally as each is presumed to fall off in the same manner.'''
# Not sure I ever finished this?  It's currently not called in the main body of the code.

#from time import asctime
#from os import linesep
import math
import numpy

import GM_InterpRoutines as IR
#import GM_logN
#import GM_N_bvals
#import GM_distance

def term3(a,b,c,alpha1,alpha2,m):
    term = alpha1*(2*a*m**(alpha2+2)/(alpha2+1) + b*m**(alpha2+1)/(alpha2))
    return term

def integrate3(a,b,c,phypars,mlo,mhi):
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    mlotemp = float(mlo)
    mhitemp = float(mhi)
    total = 0.
    for i in range(0, len(masslimits)):
        mmax = masslimits[i][1]
        mmin = masslimits[i][0]
        alpha1 = float(massmults[i])
        alpha2 = float(massindices[i])
        if (mlo > mmin) and (mhi < mmax):
            termhi = term3(a,b,c,alpha1,alpha2,mhi)
            termlo = term3(a,b,c,alpha1,alpha2,mlo)
            thisterm = termhi - termlo
            total = total + thisterm
        if (mlo > mmin) and (mhi > mmax):
            mhi = float(mmax)
            termhi = term3(a,b,c,alpha1,alpha2,mhi)
            termlo = term3(a,b,c,alpha1,alpha2,mlo)
            thisterm = termhi - termlo
            total = total + thisterm
            mhi = float(mhitemp)
        if (mlo < mmin) and (mhi < mmax):
            mlo = float(mmin)
            termhi = term3(a,b,c,alpha1,alpha2,mhi)
            termlo = term3(a,b,c,alpha1,alpha2,mlo)
            thisterm = termhi - termlo
            total = total + thisterm
            mlo = float(mlotemp)
    return total


def term4(d,e,f,g,alpha1,alpha2,m):
    term = alpha1*(3*d*m**(alpha2+3)/(alpha2+2) + 2*e*m**(alpha2+2)/(alpha2+1) + f*m**(alpha2+1)/(alpha2))
    return term

def integrate4(d,e,f,g,phypars,mlo,mhi):
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    mlotemp = float(mlo)
    mhitemp = float(mhi)
    total = 0.
    for i in range(0, len(masslimits)):
        mmax = masslimits[i][1]
        mmin = masslimits[i][0]
        alpha1 = float(massmults[i])
        alpha2 = float(massindices[i])
        if (mlo > mmin) and (mhi < mmax):
            termhi = term4(d,e,f,g,alpha1,alpha2,mhi)
            termlo = term4(d,e,f,g,alpha1,alpha2,mlo)
            thisterm = termhi - termlo
            total = total + thisterm
        if (mlo > mmin) and (mhi > mmax):
            mhi = float(mmax)
            termhi = term4(d,e,f,g,alpha1,alpha2,mhi)
            termlo = term4(d,e,f,g,alpha1,alpha2,mlo)
            thisterm = termhi - termlo
            total = total + thisterm
            mhi = float(mhitemp)
        if (mlo < mmin) and (mhi < mmax):
            mlo = float(mmin)
            termhi = term4(d,e,f,g,alpha1,alpha2,mhi)
            termlo = term4(d,e,f,g,alpha1,alpha2,mlo)
            thisterm = termhi - termlo
            total = total + thisterm
            mlo = float(mlotemp)
    return total


def cond1(xtuple,i,xval):
    if (xtuple[i] < xval) and (xtuple[i+1] > xval):
        return True
    else:
        return False

def cond2(xtuple,i,xval):
    if (xtuple[i+1] < xval) and (xtuple[i] > xval):
        return True
    else:
        return False

def borderlands(i,length):
    ends = False
    if i == 0:
        i1 = i
        i2 = i+1
        i3 = i+2
        ends = True
        ilist = (i1,i2,i3)
    if i == length-2:
        i1 = i-1
        i2 = i
        i3 = i+1
        ends = True
        ilist = (i1,i2,i3)
    if (i > 0) and (i < length-2):
        i1 = i-1
        i2 = i
        i3 = i+1
        i4 = i+2
        ilist = (i1,i2,i3,i4)
    return (ends, ilist)
    
def xytuple(xval,xtuple,ytuple,flip=False):
    # Set flip to True if the sense of xtuple is inverted
    if flip == False:
        cond = cond1
    if flip == True:
        cond = cond2
    length = len(xtuple)
    ends = False
    for i in range(0, length-1):
        test = cond(xtuple,i,xval)
        if test == True:
            (ends, ilist) = borderlands(i,length)
            (i1,i2,i3) = (ilist[0],ilist[1],ilist[2])
            x1 = xtuple[i1]
            y1 = ytuple[i1]
            x2 = xtuple[i2]
            y2 = ytuple[i2]
            x3 = xtuple[i3]
            y3 = ytuple[i3]
            if ends == True:
                (a,b,c) = IR.quadratic(x1,y1,x2,y2,x3,y3)
                yval = a*xval**2 + b*xval + c
                valtuple = (a,b,c)
            if ends == False:
                i4 = ilist[3]
                x4 = xtuple[i4]
                y4 = ytuple[i4]
                (d,e,f,g) = IR.cubic(x1,y1,x2,y2,x3,y3,x4,y4)
                yval = d*xval**3 + e*xval**2 + f*xval + g
                valtuple = (d,e,f,g)
    return (yval, valtuple)

def integrate(phypars,mlo,mhi,j,nMvs):
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    (ends, jlist) = borderlands(j,nMvs)
    (j1,j2,j3) = (jlist[0],jlist[1],jlist[2])
    m1 = MvMass[j1,1]
    M1 = MvMass[j1,0]
    m2 = MvMass[j2,1]
    M2 = MvMass[j2,0]
    m3 = MvMass[j3,1]
    M3 = MvMass[j3,0]
    mmid = (mlo+mhi)/2.
    if ends == True:
        (a,b,c) = IR.quadratic(m1,M1,m2,M2,m3,M3)
        slope = 2*a*mmid + b
        thisterm = integrate3(a,b,c,phypars,mlo,mhi)
    if ends == False:
        j4 = jlist[3]
        m4 = MvMass[j4,1]
        M4 = MvMass[j4,0]
        (d,e,f,g) = IR.cubic(m1,M1,m2,M2,m3,M3,m4,M4)
        slope = 3*d*mmid**2 + 2*e*mmid + f
        thisterm = integrate4(d,e,f,g,phypars,mlo,mhi)
    return (thisterm, slope)


def comparable(simpars, phypars):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    # Note: I use 'L'uminosity for many things that would normally be 'M'agnitude throughout this function b/c otherwise I would continually confuse them with 'm'ass
    nLs = len(MvList)
    nMvs = len(MvMass)
    # Generate comparable values from Binney&Merrifield and Kroupa distributions
    Mv_comp_vals = []
    BMlist = []
    Klist = []
    for i in range(0, nLs):
        # 1st check that you are in the range where comparison can be done
        L = MvList[i]
        nL = LumFuncList[i]
        if (i == 0) or (i == nLs-1):
            continue
        # Generate comparable values from Binney&Merrifield distribution
        L1 = MvList[i-1]
        nL1 = LumFuncList[i-1]
        L2 = MvList[i]
        nL2 = LumFuncList[i]
        L3 = MvList[i+1]
        nL3 = LumFuncList[i+1]
        (a,b,c) = IR.quadratic(L1,nL1,L2,nL2,L3,nL3)
        # Limits of integration
        Llo = (L1 + L2)/2.
        Lhi = (L2 + L3)/2.
        # 2nd check that you are in the range where comparison can be done
        if (Lhi >= MvMass[0,0]) or (Llo <= MvMass[0,-1]):
            continue
        Mv_comp_vals.append(L)
        # Integrate, record
        BMval = a*Lhi**3/3. + b*Lhi**2/2. + c*Lhi - (a*Llo**3/3. + b*Llo**2/2. + c*Llo)
        BMlist.append(BMval)
        # Generate comparable values from Kroupa distribution
        # Start by finding lower/upper m corresponding to lower/upper L
        # NOTE that Llo -> mhi and vice versa; a sign is flipped later to compensate
        (mlo, valtuplelo) = xytuple(Lhi, MvMass[:,0], MvMass[:,1], flip=True)
        (mhi, valtuplehi) = xytuple(Llo, MvMass[:,0], MvMass[:,1], flip=True)
        mlotemp = float(mlo)
        mhitemp = float(mhi)
        # Now integrate
        totalKval = 0.
        for j in range(0, nMvs-1):
            mdn = MvMass[j,1]
            mup = MvMass[j+1,1]
            thisterm = 0.
            if (mlo > mdn) and (mhi < mup):
                option = 1
                delm = mhi-mlo
                (thisterm,slope) = integrate(phypars,mlo,mhi,j,nMvs)
            if (mlo < mdn) and (mhi > mup):
                mlo = float(mdn)
                mhi = float(mup)
                option = 2
                delm = mhi-mlo
                (thisterm,slope) = integrate(phypars,mlo,mhi,j,nMvs)
                mlo = float(mlotemp)
                mhi = float(mhitemp)
            if (mlo > mdn) and (mhi > mup) and (mup > mlo):
                mhi = float(mup)
                option = 3
                delm = mhi-mlo
                (thisterm,slope) = integrate(phypars,mlo,mhi,j,nMvs)
                mhi = float(mhitemp)
            if (mlo < mdn) and (mhi < mup) and (mhi > mdn):
                mlo = float(mdn)
                option = 4
                delm = mhi-mlo
                (thisterm,slope) = integrate(phypars,mlo,mhi,j,nMvs)
                mlo = float(mlotemp)
            if thisterm > 0:
                print 'Llo', Llo, 'Lhi', Lhi, 'mlo', mlo, 'mhi', mhi
                print 'An integration element was unexpectedly positive. thisterm', thisterm, 'slope', slope, 'delm', delm, 'option', option
            # Here's the sign flip!
            totalKval = totalKval - thisterm
        Klist.append(totalKval)
    return (Mv_comp_vals, BMlist, Klist)


def fit(BMlist, Klist):
    BMlist = numpy.array(BMlist)
    Klist = numpy.array(Klist)
    KKlist = Klist**2
    KBMlist = BMlist*Klist
    KK = KKlist.sum()
    KBM = KBMlist.sum()
    c = KBM/KK
    return c
